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Abstract —Time Difference of Arrival (TDOA) is widely used 
in wireless localization systems. Among the enormous approaches 
of TDOA, high resolution TDOA algorithms have drawn much 
attention for its ability to resolve closely spaced signal delays in 
multipath environment. However, the state-of-art high resolution 
TDOA algorithms still have performance weakness on resolving 
time delays in a wireless channel with dense multipath effect. 
In this paper, we propose a novel TDOA algorithm with super 
resolution based on a multi-dimensional cross-correlation func¬ 
tion: the Volume Cross-Correlation Function (VCC). Theoretical 
analyses are provided to justify the feasibility of The proposed 
TDOA algorithm, and numerical simulations also show an 
excellent time resolution capability of the algorithm in multipath 
environment. 

Index Terms —Time Difference of Arrival, Volume Cross- 
Correlation Function, super resolution, multipath environment 


I. Introduction 

Among the tremendous amount of source localization tech¬ 
niques Ill2l3l . TDOA based techniques are widely used in 
wireless communication im indoor microphone positioning 
a, wireless sensor network 0, passive localization system 
l7l8l . and sonar J9|. Since traditional TDOA methods, such 
as the Generalized Cross-Correlation algorithm (GCC) ED, 
have limited time resolution and can not resolve the TDOA 
of multipath signals with close delays, many high resolution 
TDOA algorithms have been proposed recently to deal with the 
scenario where signals from different paths have close delays. 

There are mainly three branches of high resolution TDOA 
algorithms: one is the optimal maximal likelihood (ML) time 
delay estimators using techniques like expectation maximiza¬ 
tion (EM) Q2, or importance sampling E21; another 
branch is the super resolution TDOA algorithms based on 
subspace methods 1141151161 ; the third branch is the high 
resolution TDOA estimation methods using sparse recovery 
algorithms based on l\ optimization 1 1171181 . Except for those 
main branches, some delay estimation techniques that have 
super resolution and ability of dealing with multipath envi¬ 
ronment, such as the technique of time delay estimation from 
low-rate samples over a union of subspaces Il9l20l can also 
be adapted to TDOA estimation. 
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As we know, improving the time resolution and enhancing 
the ability of identifying each multipath TDOA are two major 
tasks concerned in design of TDOA techniques. In this paper, 
we are going to propose a highly efficient TDOA algorithm, 
which has strong ability to resolve multipath TDOAs, based on 
a novel multi-dimensional cross-correlation function, named 
the Volume Cross-Correlation function (VCC). This VCC 
function takes two matrices (which represent subspaces), in¬ 
stead of two vectors, as arguments. It calculates the geometri¬ 
cal volume of the high dimensional parallelotope spanned by 
column vectors of these two matrices. It can be regarded as a 
generalized distance measure of the two subspaces spanned by 
columns of each input matrix. In our method, the received sig¬ 
nal is formulated as deterministic signal with unknown linear 
subspace structure contaminated by random noise. Then this 
unknown subspace is extracted from noise through singular 
value decomposition of some data matrix, such procedure is 
actually a denoising process commonly seen in modern signal 
processing. Afterwards the VCC function is calculated with 
inputs being the basis for the estimated subspace. Finally the 
corresponding TDOA estimation is indicated by the zeros (or 
equivalently, the peaks of its reciprocal) of the VCC function. 

In order to analyze the performance of the proposed TDOA 
algorithm, we choose the passive localization system as a 
typical application scenario. In our analysis, the received 
signals commonly encountered in passive localization systems 
are divided into two different categories: the slowly changing 
subspace signal and the fast changing subspace signal. The 
slowly changing subspace signal means the subspace structure 
of the signal remains unchanged during the time interval of 
a large amount of observations. As for the fast changing 
subspace signal, contrast to the term ’’slowly changing”, it 
refers to the circumstance that the subspace structure are 
changing among different observations; therefore there is only 
a single observation available to estimate the current signal 
subspace. The two signal categories will cover most wireless 
signals encountered in passive localization systems. 

The rest of this paper will be as follows. In section II, we 
give the problem formulation, as well as the definition of VCC 
function. In section III and section IV, we propose and analyze 
our proposed TDOA algorithm based on two categories of 
signals, respectively. The performance of our TDOA method 
is demonstrated through numerical simulations in section V. 

II. Preliminary Introduction 
A. TDOA estimation in multipath environment 

In a typical TDOA-based localization system, due to the 
complicated environment where buildings and vehicles may 
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lead to significant scattering of wireless signals, there might 
be dense multipath effect in the wireless channel. The received 
multipath signals will be: 

x i{t) = ^2 ai,hs(t - n,^) + wi(t), (1) 

h=i 
1/2 

X 2 (t) = ^2 0 t 2 ,l 2 s(t ~ T 2 ,l 2 ) + W 2 (t), ( 2 ) 

h=l 


The Volume Cross-Correlation (VCC) function of two given 
matrices X\ £ <C nxdl and X 2 £ C nxd2 is defined as 


vcc(X 1; X 2 ) 


voldi+d 2 ([^i ; ^ 2 ]) 

vold^XOvoU^)’ 


(5) 


where [Xi,X 2 ] means putting the columns of matrices X-\ 
and X 2 together into one matrix, and vold(X) denotes the 
geometrical volume of matrix X £ <C nxd with dimension d 
(d < n). It is defined as lf24l : 


where a\^ 1 and o: 2 / 2 are the propagation gains (also known 
as the channel coefficients) of the Zith (or l 2 ) path along 
which the signal transmitting from source to receiver 1 and 2, 
respectively, and t 2 j 2 represents the corresponding path 
delays, L \ and L 2 are the number of channel paths. w±(t) and 
w 2 (t) are noises. 

As we know, the task of TDOA estimation is to determine 
the difference of time delays of the received signal from 
different receivers. However, from 0 and 0 , it can be seen 
that in a multipath channel, there are theoretically multiple 
TDOAs which can be resolved. Denote these TDOAs by 

Axil ,Z 2 • T 2 Z 2 Tt,ii ; ^1 I) 1 d/1, l 2 I) T dj 2 , (3) 

we call these multiple TDOAs as multipath TDOA. Although 
in source localization systems, the direct path TDOA is the 
only concerned, which is At^i = r 2 ,1 — Tip, precise estima¬ 
tion of the direct path TDOA Ar-i.i actually requires resolution 
of every multipath TDOA in 0 . In other words, because the 
channel path delays and propagation gains are basically un¬ 
known at the receivers, we cannot tell the difference between 
direct path TDOA and other indirect path TDOAs merely 
from the received signals. Therefore, we need to resolve every 
mutipath TDOA, before we pick the direct path TDOA and 
continue the localization process. From this point of view, the 
primary goal of TDOA localization in multipath environment 
is to precisely resolve every multipath TDOA shown in (|3j. 

B. The Volume Cross-Correlation Function 

The basic relationship between linear subspaces are gener¬ 
ally described by principal angles Ell. The principal angle is 
defined as: 

Definition 1: Consider linear subspaces Aj and X 2 , with 
dimensions dim(A’ 1 ) = di,dim(A , 2 ) = d 2 , denote m = 
min(di, d 2 ). The principal angles between subspaces Xi and 
X 2 , denoted by 0 < 9\ < • • • < 9 m < 7r/2, are defined 
recursively as 

COS 0.; = max ujVi, 

u i ex 1 ,v i &x 2 

subject to ||«i|| 2 = ||t?i|| 2 = 1, (4) 

ujuj = 0 , vj Vj = 0 , 

where i = 1, • ■ ■ m, j = 1,• ■ • ,i — 1. 

The principal angle is an important mathematical tool to 
depict the relationship between subspaces. Except for playing 
a key role in deriving the geodesic distance ll22l for Grassmann 
manifold [22] [23), principal angle are also used to define var¬ 
ious distance metrics of linear subspaces []23l . The proposed 
VCC function in this paper is related with the principal angle. 


d 

vold(-X’) := Oi, (6) 

i—l 


where o\ > a 2 > ■ ■ ■ > crd > 0 are singular values of 
matrix X. Indeed, vold(X) is the geometrical volume of d 
dimensional parallelotope spanned by the column vectors of 
matrix X. 

The relation between volume and principal angles is de¬ 
scribed by the next proposition from ED: 

Proposition 1: Consider two linear subspaces X\ and X 2 
in Ht N , their dimensions are dim(A’i) = di,dim(A 2 ) = d 2 , 
and their basis matrices are X\ £ R Nxdl and X 2 £ R Nxd2 , 
then we have 


VOlrfi+rf 2 ([-^li -^2]) 

voldj (2fi) vold 2 (X 2 ) 


min(<ii ,^2) 

l[ sin 6j{X 1 , X 2 ), (7) 

2=1 


where 0 < 9j(Xi,X 2 ) < 2tr,l < j < min(di,d 2 ) are the 
principal angles between X\ and X 2 . 

0 just indicates that the VCC function 0 is actually 
the product of sines of the principal angles between sub¬ 
spaces X\ and X 2 . Therefore 0 can be regarded as a kind 
of distance measure of subspaces X\ and X 2 . Intuitively, 
if the subspace A) and X> are linearly dependent, then 
dim(Ai fj X 2 ) > 0. According to the definition of principal 
angles, there must be a vanishing principal angle 9j(X i, X 2 ). 
That is vold 1+ d 2 ([JVi,X 2 ])/vol dl (Xi) vol d2 (X 2 ) = 0. 

On the other hand, if X\ is perpendicular to X>, then 
vold 1+ d 2 ([Xi, X 2 ])/voldjCXi) vold 2 (X 2 ) = 1 holds obvi¬ 
ously. As a matter of fact, VCC function measures the extent 
of linear dependency between subspaces ||25 1 . and will be used 
to derive our TDOA algorithm. 


III. Estimating the TDOA of slowly changing 
SUBSPACE SIGNALS USING VCC FUNCTION 
A. The slowly-changing subspace signal 

In the passive localization system, information about the 
wireless channel as well as the source signals are generally 
unknown by the receivers. Hence TDOA technique is quite 
suitable for this kind of localization system B26I27B . Firstly, we 
focus on the category of signals that have a slowly changing 
subspace structure. 

A typical type of signals we encounter in passive local¬ 
ization systems are radar signals radiated by non-cooperative 
radar transmitters. The common pulse radar waveform can be 
expressed as the following expression: 

+oo 

s(t)= ^2 '/Psg{t-mTp), (8) 
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where P s > 0 is the transmitting power of the radar, and 
g(t ) £ C is the general form of the radar pulse waveform, T p 
is the pulse repetition interval (PRI). 

We consider the received signal in a sigle PRI, then in 
the multipath environment, the received signal from the i’th 
receiver ( i = 1, 2, • • •) would be 
Li 

Xj(t) = y saiMy/Psdjt-Tiji) £e(0,Tp), (9) 

li =1 

where and T, j t are multipath channel coefficients and 
path delays, vj, (t) is the Gaussian white noise of Pth receiver. 
After sampling the received signal with the rate 1/T S , then (0 
can be written into 

Li 

yi{kT s ) = ^ a i ,i i '/P~ s g(kT s - t*,j 4 ) + Wi(kT s ). (10) 
h=1 

Rewrite (HI into a vector form, we have 

Vi — \/PsG/OLi Wij i — 1, 2, • • • , (11) 

where := [y»(0), yi{T s ), ■ ■ ■ yi((N - 1)T S )] G C N , and 

G, = [<7»,i, <7i, 2 , • • ■ ,9i,Li] S C Wxi % (12) 

where g i t := [c/(0 -T s -r^J, ff(l -T s -t^J, ■ • • ,g((7V-l)- 
2 n s-A,i i )] T ,( = l,-" ,Li- The vector a * = ■ ■ ■ a itLi ] G 

C Li is the channel coefficient vector composed of the malti- 
path channel gains, and Wi is the noise vector. As is shown in 
CD- the received radar signal in a multipath channel generally 
has a deterministic subspace structure with the corresponding 
subspace span(G), i.e., spanned by different time-shift ver¬ 
sions of radar waveform g(t). 

Except for radar signals, the common linearly modulated 
wireless communication signals such as DS-CDMA, OFDM, 
QAM, and others that carry symbols on some periodic pulse 
shapes, can also be modeled as the signal with a subspace 
structure in O 1281291301 . The subspace signal structure 
is mainly related with the channel’s path delays As a 
matter of fact, channel delays are caused by different distances 
between receivers and signal sources (or reflective objects), 
and generally signal sources and reflective objects don’t have 
extremely high velocities, therefore channel delays can be 
generally seen to be constant in a short time. On the other 
hand, the wireless channel’s path gains a fluctuates with 
time, which is caused by channel fading effect. This fact 
also means that for a time interval long enough for the 
receiver to obtain relatively large samples of the received 
signal (according to chapter 2 of ED, the time scale of this 
interval can be up to 20s in a typical channel scenario), these 
sample data can be formulated as 

y^=G ia ^+w^, j = 1,2, ■ ■ ■ , (13) 

where j indicates different observations at different time, or 
snapshots. Although the channel coefficient vector a.G) i s 
fluctuating with different j, the subspace structure determined 
by matrix G,, will remain almost unchanged. We call this cat¬ 
egory of signals the slowly changing subspace signal, meaning 
that the subspace structure in (ITD changes slowly with time, 
and can be treated as invariant during the observation interval. 


The subspace structure span(Gj) is unknown to the re¬ 
ceivers, but can be estimated from multiple observation data 
like (HD. Because a typical radar transmits a pulse waveform 


Radar 



Fig. 1: Getting multiple observations of radar signals 


repeatedly with a PRI of T p , we can receive multiple snapshots 
of these multiple radar pulses as in (fl3l > according to T p , which 
can be estimated by various PRI identification techniques 
H32I331I , Afterwards the corresponding signal subspace is esti¬ 
mated using the well-known subspace methods like MUSIC, 
ESPRIT, etc. The process of obtaining multiple observations 
and estimating subspace span(Gj) can be demonstrated in 
figure Q] the gradient change of the background color in figure 
(D represents the fluctuation of channel coefficients a they 
are reasonably assumed to take independent values among 
different pulses’ durations. 

Denote these multiple data of the received radar signal by 


V 


(i) 


■V 


(m) 


(14) 


they are then used to evaluate the sampled covariance matrix 

1 m 

= ( 15 ) 

TIL 

3 = 1 

According to the well-known subspace methods, the signal 
subspace can be estimated through eigen-decomposition of Ri. 
We can firstly evaluate the eigenvalues and eigenvectors of Ri, 
then we estimate the dimension of the signal subspace (if the 
matrix Gi is full rank, the dimension will be /. Jj]) by analyzing 
the distribution of eigenvalues; and finally we can separate the 
eigenvectors of Ri into bases for signal subspace as well as 
bases for noise subspace. The bases for signal subspace are 
eigenvectors with respect to the L, largest eigenvalues. As a 
result, we can write Ri into 


Ri = Ui^Ai^U^g 


Ui. 


A U H 


(16) 


where the matrix C/) jS is the estimated basis matrix for the 
signal subspace span(Gi). It has been proved that span([/i :S ) 
approximates the signal subspace span(Gi) asymptotically for 
a sufficiently large m I34l35l . Then we will use the estimated 
basis of signal subspace, i.e., Ui tS , to estimate TDOA using 
VCC function. 


B. TDOA estimation using VCC function for slowly changing 
subspace signals: Algorithm 

Algorithm 1. 

'in following analysis, we will assume this full-rankness to be satisfied. 
Indeed, a sufficient condition will be given to ensure this assumption. 
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For At G ((-N + 1 )T S , (N - 1 )T S ), 

1) Obtain delayed multiple observation data (for 

j = I,-" ,m) 

„,U)A At ] _ 
y i — 

[wi O) (0 • T, - At), ■ • • , y[ j) ((N - 1 )T S - Ax)]' 
and non-delayed observation data 
vi j) = ,y a 0) ((JV-i)T.)] T . 

2) Calculate the sampled covariance matrices 
r[ A t] and R 2 according to (fl5l ). 

3) Perform eigenvalue decomposition of the sam¬ 
pled covariance matrix R Ar (and R/)', 

4) Estimate the dimensions of signal subspaces L\ 
(and L 2 ) based knowledge of the eigenvalues; 

5) Extract the basis for the signal subspaces 
U[ A S T] G C NxLl (and U 2 , s G C NxL2 ) from 
the eigenvectors corresponding to the L\ (and 
L 2 ) largest eigenvalues of r) A t ^ (and R 2 ): 

6 ) Calculate the reciprocal of VCC function: 

r»i(Ar) := 1/ vol Ll +L 2 ([U[ A a T] , U 2 , s ])- 


Find the peaks of r ml (Ar) and the corresponding At. 

Below are some supplementary remarks: 

Remark 1. In our algorithm, the dimensions of signal 
subspaces are actually the number of channel propagation 
paths, i.e., Li in ©, which are generally unknown at the 
receiver. Therefore the dimensions have to be estimated first. 
Luckily there are various methods to estimate the dimension 
of signal subspace from a sampled covariance matrix, such as 
Akaike Information Criterion (AIC) |36) Minimum Descrip¬ 
tion Length (MDL) |37| , Bayesian Information Criterion (BIC) 
lf38l . Predictive Description Length (PDL) |f39l and so on. In 
this paper we just assume the number of channel paths Li is 
precisely estimated. 

Remark 2. The basis matrices u[ A t ^ and LL.. S are the 
estimated bases for the delayed signal subspace span(G^ Ar ^) 
and non-delayed subspace span(C? 2 ), respectively. From the 
previous analysis, we know that 


G [At] = 


[Ar 

01,1 


o [Ar 


G C 


NxL 1 


(17) 


where g [A ^ := [g (0 • T s - At - r lth ), g( 1 ■ T s - At - Tin), 
■ ■ ■ ,g((N - 1 )T S - At - Tin)] T , and 

G 2 = [g 2 ,t,--- ,g 2 , L2 ] GC“ 2 , (18) 


where g 2 ,z 2 := [ 9(0 • T s - r 2 | j 2 ), 9(1 • T s - t 2 ,i 2 ), ••• ,g((N- 
l)T s -t 2j ; 2 )] t . So the VCC function vol Li+ l 2 ([u[ At \u 2 ^ s ]) 
is approximately measuring the linear dependence of sub¬ 
spaces span(G[ Ar ^) and span(G 2 ). It is obvious that 
when At = t 2 ,; 2 — Tij lt for any pair of l\ and l 2 , 
span(G 2 Ar ^) and span(G 2 ) are linearly dependent, which 
means 1 / yoXl^+l^ U 2 , s }) will tend to infinity; while 

on the other hand, when At 7 ^ t 2 j 2 — Ti t u span(G[ Ar ^) 
and span(G 2 ) may be linearly independent, meaning that 


1/vo1l 1 + l 2 ([g| A t ^ U 2jS ]) would have a finite value. From 
this observation we could expect that r vo i(Ar) would reach 
its peak at every At = t 2 j 2 — Tp/j. We will give a theoretical 
guarantee in the next section to validate this observation. 

Remark 3. It should be noted that the reciprocal of VCC 
function r vo i (At) is actually a continuous function of variable 
At. So the final step of the proposed algorithm actually 
involves a continuous one-dimensional parameter search. In 
practice, when we are manually delaying the signal y^\ At 
can only be integer multiple of T s , i.e., At can only take 
discrete values related with T s . However fortunately, as long 
as the signal is sampled at Nyquist rate, we can get the signal 
yp)'I Ar l for any At without loss by interpolating the original 
signal until At is on the sampling grid. This interpolation can 
be a pre-process of the algorithm and is verified by simulation. 

Remark 4. Theoretically, we are expecting to resolve every 
multipath TDOA, i.e., 

At / 2 li . T 2 i 2 T[ , 

for all (1 = 1, • ■ • , Li, and l 2 = 1, • ■ • , L 2 . It is another inde¬ 
pendent problem about how can find the TDOA corresponding 
to the direct channel path, i.e., identifying Ati.i = T 2 ,i — ti,i 
among all of the multipath TDOA, which is known as the dis¬ 
ambiguation of TDOA 1401411 . This topic won’t be discussed 
in this paper and will be left for a future work. 


C. A Theoretical Guarantee for Algorithm 1 

Denote the auto-correlation function of radar waveform q(k- 

T s ) by 

N-l 

R a( r ) '■= 9(kTs) ■ g*(kT 3 - t), 

k =0 


The auto-correlation function is well known as the ambiguity 
function of a radar waveform along the zero-Doppler axis, it 
is an important characteristic of the radar pulse waveform. 
We will show here that, the performance of the proposed 
algorithm is theoretically guaranteed, and related with the 
auto-correlation function of the transmitted radar waveform. 
A lemma is need firstly: 

Lemma 1: Consider the matrices G 7 and G 2 in (fT71 > and 

m, if 


|77 g (r)| 1 

\R g (0)\ L 1 +L 2 -l 


for |t| > At*, 


then rank(Gj Ar ^) = Li,rank(G 2 ) = L 2 . Here 


(19) 


At* := min{ATi, AT 2 }, (20) 

and 


An := min Inn - n ; ' 

At 2 := min |t 2 ,i 2 - r 2 i 2 /|. 

l<h^h'<L 2 

Lemma Q] provides a sufficient condition for the full¬ 
rankness of matrices G 7 and G 2 . Under the condition ( I I b| ) 
of Lemma[I] all the multipath TDOAs t 2 ,z 2 — Tin will theoret¬ 
ically be resolvable by our algorithm, the following theorems 
will describe the theoretical behavior of your algorithm when 
Lemma Q] holds: 
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Theorem 1: When Lemma Q] holds, consider the matrices 
G^ A ^ and G 2 in ( ITTl i and ([T 8 K when for any 1 < h < L\ 
and 1 < I 2 < L 2 , At 7 ^ T 2 ,z 2 — ti^, given any parameter /i 


satisfying 

L 1 +L 2 -I’ 


(21) 

if 


km sfor h - aTm,n ’ 

(22) 

then we have 


vol Ll+L2 ([C/' Arl ,G 2 ]) > (1 — e) i/2 , 

(23) 

where 


Ar min := min |Ar - t 2 ,z 2 + 

i<Zi<Li,i<z 2 <l 2 

(24) 

L = min {L \. L 2 } and the parameter e is 


L 1 L 2 ■ / 1 2 

[1 - (L, - 1M1 - (£2 - 1)4 

(25) 

Theorem 2: If there exist If and If, such that At 
n,!*, then 

= - 

vol Ll+i2 ([[/j AT] ,[/ 2 ])=0. 

(26) 


The matrices u[ A t ^ and U 2 here are orthogonal basis matrices 
for subspaces span(G^ AT ^) and span(G 2 ). 

Theorems Q] and [2] provide sufficient conditions for the 
feasibility of our proposed TDOA algorithm. Theorem[2] shows 
that, when there exist If and If, such that 

At = t 2 ,;* - ti,j., 

the VCC function will certainly be zero. This means 
1/ vo 1 Li+L2 ([[/| At] , C/ 2 ]) reaches infinity when At equals the 
TDOA corresponding to path delays T 2 j* and t\ j -. On the 
other hand, according to theorem |T] when 

At 7 ^ t 2 ,i 2 - t 1m , 

as long as the condition (l 22 l> holds, the volume function 
volL 1 -t-L 2 ([C/| Ar ^, C/ 2 ]) will be non-zero (because e is surely 
less than 1 for 0 < p < l/(Li + L 2 — 1 )) as shown in 
(|23K Because Theorems Q] and [2] ensure the different values of 
vo1l 1 _)_l 2 ([C/| At ^ C/ 2 ]) when At is under different conditions, 
thus guarantee the feasibility of our proposed algorithm for 
identifying multipath TDOA. 

It is worth noticing that we are using the reciprocal of 
VCC function, i.e., 1/ vo1l 1 +l 2 ( [C7’| At ^ , C/ 2 ]) to identify the 
multipath TDOA, because it approximately reaches infinity 
when At equals one TDOA, and otherwise, remains a finite 
value. This procedure is similar to the reciprocal technique 
commonly known in MUSIC method, thus we could conjec¬ 
ture a similar super-resolution characteristic when identifying 
TDOA. Numerical simulations in section V also show that the 
VCC function does reveal sharp peaks at the corresponding 
TDOA locations, thus has super-resolution. 

If we take a close look at (fl9l > and d22k it is obvious that the 
behavior of our algorithm is related with the autocorrelation 
function of the radar waveform g(kT s ), i.e., R 9 {t). Basically, 


the capability of identifying different TDOAs is determined by 
the autocorrelation function R g (r). Typically |i? g (T)| reaches 
its maxima at t = 0; and \R, r/ (T)\ will eventually (or oscilla- 
torily) decrease as |t| increases. So if the radar waveform’s 
autocorrelation function has just one sharp and narrow peak at 
t = 0 , the proposed algorithm will be able to identify different 
TDOAs that are close to each other; while on the other hand, 
when | R g ( t )| drops slowly as |t| increase, then for example, 
if there are If and If' such that — T lt ;*/ are too small 
to satisfy < 0 , TDOAs related with ti ;* and t \;»/ might no 
longer be identifiable. 

As we have discussed, signals with different autocorrelation 
functions will cause different TDOA estimation performance 
of our algorithm. Because the width and sharpness of mainlobe 
of auto-correlation function is actually related with the corre¬ 
sponding signal’s bandwidth, wideband radar signals would 
bring better precision in our TDOA estimation algorithm, 
simulation will be carried out to support this conclusion. From 
this point of view. Theorems [ 2 ] and Q] coincides with traditional 
theoretical analyses of TDOA methods like GCC, where the 
bandwidth of signals is always a key factor influencing the 
TDOA esimation precision. 


IV. Estimating the TDOA of fast changing 

SUBSPACE SIGNALS USING VCC FUNCTION 


A. The fast changing subspace signal 

Contrast to the slowly changing subspace signal model, 
there are also a large category of signals that don’t have a 
steady subspace structure as in (fill . For example, in passive 
localization systems, FM radio transmitters, TV broadcast 
stations are usually the signal sources to localize, or are used as 
the illuminators-of-opportunity to localize a reflective target. 
Because this category of signals are randomly varying with 
time and have no repeating waveforms, we cannot get multiple 
observations of the received signal as in (IT3l > that have the 
same subspace structure. This category of signals is called 
the fast changing subspace signal. In this case, the received 
baseband signal from the / 1 h receiver is in the form of: 

L, 

x i(t) = ^ OLiM s {t - + Wi(t ), i = 1, 2, • • • , (27) 

h =1 


where and T* i ; i are channel’s path gain and path delay, 
respectively. The original transmitted signal s(t) can be FM, 
PSK or AM signals, etc. 

Although we cannot estimate the signal subspace from 
multiple observations as in the previous section, there is still a 
way to extract a time-dependent signal subspace from a single 
observation data. Suppose a sample rate T s , for a single obser¬ 
vation data x = [cc(0 • T s ), a;(l • T s ), ■ ■ • , x((N — 1)T s )J t , we 
can construct its Hankel matrix (also referred to as trajectory 
matrix), which is 


X = 


x{0-T a ) x(l ■ T s ) 

x(l-T a ) x{2 -T s ) 


x((K - 1 )T S ) 
x(KT a ) 


(28) 


Li((M- 1)T s ) x{MTs ) 


x((N-l)T s ) J 


where 1 < M < N,K = N—M+ 1. The left singular vectors 
of the Hankel matrix X are known contain important informa¬ 
tion about the signal x l42l . Therefore the subspace spanned 
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by a subset of these left singular vectors is called ’’the signal 
subspace” (generally the left singular vectors corresponding 
to larger singular values would be chosen). As a matter of 
fact, this signal subspace extracted from the Hankel matrix 
has been used to perform noise reduction, signal forecasting, 
and change point detection, etc I43l44l . The Hankel matrix 
technique can be used to analyze a wide variety of signals, 
like wireless signals, seismologic, meteorological, geophysical 
time series as well as economic time series. Because no 
statistical assumption concerning the signal is needed while 
performing the subspace extraction from Hankel matrices, 
this methodology is suitable to deal with the fast changing 
subspace signal and develop our TDOA algorithm. 

In our setting, at the receiver i (i = 1,2), given 

the sampling rate 1 /T a , the sampled signal vector is 
[xi(0), x(T s ), ■ ■ ■ ,Xi((N — 1) • T S )] T with length N, denote 
the corresponding Hankel matrix as in (1281) by X, £ C M/K , 
Xi can be further written as 


x l ^J2 a ^ slTi ' h] + w l , 

li = 1 

where 1 < M < N. K = N — M + 1, and 
:= 

s(0 • T s - ••• s((K - 1)T S - Ti^) 

s(l-T s — t m .) s(KT s -T i: i.) 

_ s((M-l)T s -T iM ) ■■■ s((N-l)T s -r Mi ) . 


(29) 


(30) 


is the Hankel matrix of the sampled transmitted signal s(kT s ) 
delayed by . Wi is the Hankel matrix of noise Wi(kT s ). 

Denote the column vectors of matrix by 

, ■ - • ,s^\ in the following analysis, we just assume 
these column vectors to have equal length, or instant power, 
i.e.. Us' 1 / || 2 = = || || 2 . This assumption is quite 

general because signals we are interested in can be regarded 
as stationary so the power is treated as time-invariant during 
the time of a single observation. The full algorithm for TDOA 
estimation of fast changing subspace signal is given in the 
following. 


B. TDOA estimation using VCC function for the slowly chang¬ 
ing subspace signals: Algorithm 


Algorithm 2. 

For At £ ((-Af + 1)T S , (M - l)T a ), 

1) Construct delayed Hankel matrix 


xi(0T 3 —At) x\((K — 1)T S — At) 

xi(l-T s — Ar) ■■■ xi(KT a — At) 


_ X!((M - 1)T s - At) X1 ((N - 1)T S - At) 

and non-delayed Hankel matrix 

x 2 = 

' X 2 (0 -Ts) ■■■ x 2 ((I<-l)T a ) ' 

x 2 (l -T a ) ■■■ x 2 (KT s ) 


x 2 ((M — 1)T S ) ••• x 2 ((N — l)T a ) 


2) Compute singular value decomposition of 
and X- 2 , then we choose a subset of their left 
singular vectors, i.e., 

u^ 1 , ■ ■ • , u [ *k[ , 1 < A'l < min(M, K) 

and 


it 2 ,i, ■ ■ • , ui,k 2 > 1 < K 2 < min(M, K) 

which correspond to the singular values crip > 

Cl 2 , ■ • • , > (7i,Ki > 0 of matrix and singular 

values 0 - 2,1 > 0 - 2 , 2 ) ■ ■ ■ , > cr 2 ,A ' 2 > 0 of matrix X 2 - 
Then the matrices 

u[ A t] := [4 A i Tl ,--,4 A x!]ec™ 


and 


U 2 := [- 1 * 2 , 1 , ■ ■ ■ - 1 * 2 , K 2 \ £ C NxK2 


are basis matrices for the signal subspaces of 
receiver 1 and 2 . 

3) Calculate the reciprocal of VCC function: 


rw(Ar) := 1/ yo\ Ki+K2 {[u[ At] ,U 2 }), 

Find the peaks of r vo i(Ar) and corresponding At. 

Remark 1. Similarly r vo i(Ar) is also a continuous function 
of variable At. And because At can not take continuous 
values in practice, a similar interpolation step can be used 
to get the signal with any delay At before we construct the 
Hankel matrix, simulation will also be provided to show the 
effect of interpolation. 

Remark 2. There are two important parameters when con¬ 
structing the Hankel matrix, i.e., the dimensions M and K. It 
is difficult to choose these two dimensions in order to meet 
different requirements in diverse applications 1441 . so in this 
paper we just choose these two dimensions empirically based 
on the experiments and simulations. 

Remark 3. The eigenvector extraction procedure can be 
regarded as both feature extraction and noise reduction. An 
important parameter affecting the extraction of signal subspace 
and calculation of VCC function is the dimension of signal 
subspaces, i.e., K t . This parameter will be also determined 
empirically. Actually, in the numerical simulation which will 
be shown in the next section, Ki is chosen to be 3 times of 

U 


C. Theoretical Guarantee of Algorithm 2 

Theorem 3: Denote the auto-correlation function of trans¬ 
mitted signal s(kT s ) by R s (t), consider the Hankel matrices 
and X-> in Algorithm 2, When for any 1 < l\ < 
Li,l < l 2 < 2 * 2 , At ^ T 2,; 2 — Ti,ii- gi yen an y parameter 
/* satisfying 

0<p<^, (31) 

< At, for |t| > AT min , (32) 

\rl s {v)\ 

then the volume function in algorithm 2 must satisfy 

v °lAi+if 2 ([l^j Ar] , ty 2 ]) > (1 - e 2 )W 2 


(33) 
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where 

Ar min := min{An, At 2 , An i2 }, (34) 

and 

An = min Inn - Dip + (k - k')T s |, 

l<h^h'<Ly 
1 <k,k‘<K 

At 2 = min |t 2 j z 2 - t 2 ^ + (k - fc')T s |, 

1 < Z 2 t "^2 < 1/2 
1 <k,k'<K 

An 2 = min |n,«i + Ar - r 2 ,; 2 + (A: - fc')T s |, 

l<Zi<Li,l<J 2 <i 2 ’ 

1 <k,k'<K 

Kmin ■= minjA'i, A" 2 }, and the parameter 

£ = C • n, (35) 


here 




c = 


&l,Ki &2,K 2 


(36) 


and s 2 .i 2 represent column lengths of Hankel matrices 
S'[ T i’ i i +AT l and respectively; and cr 2i #; 2 are the 

A'i ’th and AT 2 ’th largest singular values of matrices xj Ar ^ 
and X 2 . 

Theorem 4: If there exists Z* and Z|, such that Ar = r 2 .;» — 
n,;*, given any parameter /i satisfying 


0 < p < D/(L — 1), (37) 


< A' for M > A< in , (38) 

then the volume function will satisfy 

vol Kl+K2 ({u[ Ad - T ‘\ U 2 }) < (1 - 7 2 ) i W 2 ! (39) 

meanwhile, the right side of (l39l > is less than the right side of 
d33|. Here 

Ar* nin := min{An, At 2 }, (40) 


the parameters 


D = 


' <7l,K 1 ^2,K 2 ■ A - (L — 1) 1 1 

(B + B*) ■ K + 4 _ 2’ 


(41) 


and 


A 


B* K 


7 = 


i + (L — l)n n,n ff 2 ,if 2 

where L = min{Li, A 2 }, and 

_ 

[At 


• ^ 


A := 


(42) 


B := ^2 Ms[V ' \ a 2,h\s2,i 2 

h=l Z 2—1 

B * ■= Y2 lai.hls^ 1 • \ a 2,i 2 \ s 2,h, 


Theorem |3] and H] provides a similar theoretical guarantee 
for the proposed algorithm 2 as Theorem Q] and |2] does. As 
can be seen, CH and <[33 ensure the different values of 


vo\k 1 +k 2 {[u[ At \u 2 ]) when Ar as well as and r 2j ; 2 
satisfy condition (f32l) or (f38l) . which is also related with the 
autocorrelation function R s (t). 

However, there is a major difference between Theorem [4] 
and Theorem[2] When there exists some I* and 1%, such that 

Ar = r 2j ;* - n,i*, 

the VCC function vo\k 1 +k 2 {[U[ At \u 2 \) won’t necessarily 
be zero, as the result in (l39t shows. As a result, the reciprocal 
1/ vol k 1 +k 2 ([U[ At ^ , U 2 ]) would show a finite peak when <f38l) 
holds. And the parameter 7 shown in (l42l >. which is mainly 
dependent on channel’s path gains 07,7 and ct 2 i z 2 , will affect 
the sharpness of their corresponding TDOA peak at location 
Ar, then influence the precision of estimation of their TDOAs. 

Similar to the previous analysis, the conditions (l32l) and 
( 1381 ) also indicate that, as long as the autocorrelation function 
R s (r) is ’’sharp”, the performance of the proposed algorithm 
can be theoretically better. Roughly speaking, when the signal 
has an extremely sharp autocorrelation function R s (t), all the 
parameters such as /x and Ar m j n , Arj7 n can take a small 
value, meaning a better TDOA resolution and sharper peaks 
when r is at TDOA positions. This also means that the 
proposed algorithm prefers signals that has a wide bandwidth. 

V. Numerical simulations 
A. TDOA of slowly changing subspace signals 

Firstly, a demonstration of the TDOA algorithm’s out¬ 
put is given by simulation in figure [2] In the simulation, 
a linear frequency modulation (LFM) waveform is chosen 
as a typical slowly changing subspace signal, which is the 
most commonly seen waveforms in radar systems. The radar 
waveform in © is generated with a sample rate 1 MHz, 
its length are 2048, and the frequency sweeps linearly from 
50kHz to 500kHz. The multipath channel are manually 
generated, and for convenience, the multipath delays are 
chosen arbitrarily to be exactly ’’integer delays”, in other 
words, = {40T S ,75T S ,200T S } and = 

{50T S , lOOTs, 185T S , 250T S }, so that the TDOA can theoret¬ 
ically recovered by testing the VCC function upon integer 
delays of Ar, i.e., Ar = (■ • • , 0,1, 2, • • ■) • T s . The multiple 
observations in the form of O are directly generated by 
Monte-Carlo method, in which the channel coefficients af' 1 = 
[al l , • • • , a f L . J with respect to different j are generated 
independently from complex Gaussian distributions in order 
to simulate the channel fading effect. In addition, the mean 
value of |aj,i| is greater than the mean value of 107 ^ |, li > 1 , 
meaning that the direct path has a greater propagation gain 
than the reflective path. The length N of each observation 
vector is 512, and totally 512 observation data are generated. 

In the simulation, we compare our proposed TDOA al¬ 
gorithm with the publicly known super resolution MUSIC- 
Type TDOA algorithm proposed by Fengxiang Ge in tm 
because both algorithms have super resolution and can make 
use of multiple observation data. Since the simulation focuses 
on demonstrating the ability of resolving multipath TDOA, 
we just assume the dimensions of signal subspaces in both 
algorithms, i.e., the number of channel paths Li, have been 










Comparison of Music-Type algorithm and VCC function for TDOA estimation, with SNR=-10 



a ^ A 1. A i. 


Ar/T, (Discrete Time) 


Fig. 2: Comparison of Ge’s MUSIC-Type algorithm and the 
VCC algorithm 


accurately estimated. The normalized TDOA estimation results 
of both algorithms are plotted in figure [2] where the signal- 
to-noise ratio SNR, defined as the power ratio of signal 
and noise, is set to be — IOgLB. According to the simulation 
setting, there should be peaks at positions where At = 
(-150, -100, -25, -15,10,25,50,60,110,145,175,210 )-T s 
in the TDOA estimation outputs. The position of these peaks 
are labeled in the figure. As shown in figure [2] when the SNR 
is low, the MUSIC-Type algorithm fails to reveal most of the 
peaks of multipath TDOA, but our VCC algorithm can still 
show clear peaks. 

Secondly, the overall performance of both algorithms are 
given in figure [3] In the simulation, 120 independent trials 
of both algorithms with arbitrary multipath delays are carried 
out for different SNR levels. Because the output of these algo¬ 
rithms have different scale, we define the mean square errors 
(MSE) of TDOA estimation to be MSE = — 

/(At)) 2 , where t(At) is the normalized output result shown 
in figure [2] of each algorithm, and /(At) is the ’’standard 
output vector” which takes value 1 when At is at these 
multipath TDOA positions and 0 otherwise. The simulation 
results in figure [3] implies that, the proposed VCC algorithm 
outperforms Ge’s Music-Type algorithm at all SNRs; besides, 
our VCC algorithm will have better performance when source 
signals have wide bandwidths. 

B. TDOA of fast changing subspace signals 

In this part of simulation, we chose a set of real-world 
frequency modulation (FM) broadcast signals as one example 
of the fast changing subspace signal, to demonstrate the TDOA 
estimation performance of our proposed method. The FM 
signals used here are baseband signals transmitted by a real 
world radio broadcast station gathered from several remote 
located radio receivers. 

1) Real world FM signal, when only one single path exists: 
In the simulation, the FM signals of a radio station are sampled 
from two separately located radio receivers, the sample rate of 
received baseband signals is 256kHz, and the length is 4096. 

We firstly increase the original sample rate of the raw signals 
by a factor of 4 before we use them for TDOA estimation, 


MSE of TDOA Estimation (VCC and Music-Type) 



MSE of TDOA Estimation (VCC, LFM of Different Bandwith) 



Fig. 3: Comparison of the VCC algorithm and MUSIC-Type 
algorithm (above) and Performance of the VCC algorithm 
using LFM signal with different frequency range (below) 


a part of the waveform in time domain and the frequency 
spectrum of these two baseband signals are plotted in figure 
[4] From the waveform of both signals, we can see that the 
corresponding discrete time TDOA is from 14 to 16. 

In the simulation, we compared our VCC algorithm with 
the traditional GCC-PHAT method, the high resolution (i-\ reg¬ 
ularization algorithm, and also the super resolution MUSIC- 
Type algorithm by Ge. In the simulation of £i regularization 
algorithm, the power spectrum of the transmitted signal is 
required to be known, while the other three algorithms don’t 
use knowledge of the power spectrum. In our algorithm, 
the parameters N, M and K, are chosen empirically to be 
N = 544, M = 512, = K 2 = 3. The normalized 

TDOA estimation results of GCC-PHAT, i\ regularization, 
MUSIC-Type as well as VCC algorithm are shown in figure 
0 It can be seen that in a channel with only a single path, 
both our proposed VCC algorithm and Ge’s MUSIC-Type 
algorithm outperforms the traditional GCC-PHAT and the l\ 
regularization algorithms; because the latter two methods give 
a much wider peak, and also reveal too many false peaks 
except for the real TDOA peak. Although our VCC method 
and MUSIC-Type algorithm have similar super resolution 
ability, the computational complexity of our method is much 
lower. 



Fig. 4: Real world FM signals from two separated receivers 
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Fig. 5: Comparison of different TDOA techniques in a single 
path channel environment 


2) Real world FM signal, the multipath channel is manually 
simulated: In this simulation, the received signals from two 
receivers are generated as the following expression: 

yi{kT s ) = ai,i s(kT s ) + oi, 2 s((k — 60 )T S ) + ai,3 s((k — 120 )T S ), 
y 2 (kT s ) = a 2 ,is((k - 25 )T a ) + o 2 , 2 s{{k - 100)T S ) 

+ai, 3 s((k — 195)T S ). 

The s(kT s ) here is a real world original FM signal mentioned 
before, which is also one among the two signals plotted in 
figure |4] The channel coefficients = 1,2, j = 1,2,3 

are also generated to simulate a Rician fading channel, among 
these coefficients the mean value of \aip\ is greater than that 
of the other coefficients. In the simulation, the parameters N, 
M and K, are also chosen empirically to be N = 896, M = 
768, Ki = K ‘2 = 9. The TDOA estimation results of GCC- 
PHAT, t\ regularization, MUSIC-Type and our method are 
shown in figure [ 6 ] 






Fig. 6 : Comparison of different TDOA techniques in a simu¬ 
lated multipath channel environment 

As is seen, both Ge’s MUSIC-Type algorithm and our 
VCC algorithm outperforms the other two methods. However, 
Ge’s MUSIC-Type method and our VCC algorithm have their 
advantages and disadvantages at different aspect. The MUSIC- 
Type method has a much sharper peak, but fails to resolve 
every multipath TDOA, and still has some false peaks, while 


our VCC method may not have such sharp peaks, but success¬ 
fully reveals every TDOA peak precisely with no false peak. 
In addition, our VCC algorithm has much lower computational 
efficiency, because both the MUSIC-Type algorithm and the t\ 
regularization algorithm contain a convex optimization step. 


C. Simulations on non-integer multipath delays 

In the previous simulations, we assume the multipath delays 
to be exactly on the sampling grid, i.e., Ti,ijT s takes integer 
value. But actually, in reality, r^ ^/Tg does not necessary 
take integer value, as a result, the corresponding TDOA At 
won’t be integer multiple of T s either. So we’re validating 
the performance of our algorithm when there is non-integer 
multipath delays. The simulation in figure [7] uses the same 
signal model as before except for some non-integer delays, 
it shows that, with a simple interpolation of the received 
signal before we construct the signal subspace, the non-integer 
TDOA delays can still be identified. In the upper figure, 
we interpolate the sampling rate of received signal Xi(kT s ) 
by 5 times, and in the lower figure, the received signal 
Xi(kT s ) is interpolated by factor 4. As is shown, the VCC 
function correctly show peaks at the right TDOA positions. 
This also suggests that interpolation can be used as a pre- 
process before the VCC algorithm to improve the precision of 
TDOA estimation. 



0 5 10 15 20 25 30 35 40 45 50 

At/T s (Discrete Delay) 



Fig. 7: TDOA estimation with non-integer multipath delays 


VI. Conclusion 

In this paper, a super resolution TDOA estimation technique 
using the Volume Cross-Correlation function is proposed. 
This technique firstly estimates the unknown signal subspace 
from the received signal, and estimate the time difference 
through the novel VCC function, which calculates the linear 
dependency of these subspaces. We analyzed the performance 
of our TDOA estimation algorithm upon two typical categories 
of signals, i.e., the slowly changing subspace signal and 
the fast changing subspace signal. Analysis and numerical 
simulations have demonstrated that our algorithm has excellent 
capability of super resolution for TDOA estimation in a 
multipath environment. 
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Appendix 


A. Proof of Lemma [ 7 ] Theorems \T\and [2} 

Several lemmas are needed first: 

Lemma 2: Given a matrix X = [xi, • • • ,xf\ £ 
C NxL ,L < N, denote its maximal column correlation (or 
coherence) by 


H := max 


then if 


I (xi,xi')\ 

l Al' 11**112 • ||*i' lb ' 

1 

H < 


L — 1 


(43) 

(44) 


the matrix X is full rank. 
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Proof of Lemma [2} This lemma is actually the famous 
Gershgorin circle theorem, the proof can be found in com¬ 
monly seen textbooks on matrix analysis, therefore we ignore 
the proof here. 

Lemma 3: Consider the matrices X\ = 

[*i,i, • • • ,*i,lJ G C Arxil and X 2 = [* 2 ,i,--- ,x 2 ,l 2 \ G 
C A y L ' 2 . L \. A 2 < N denote the maximum column correlation 
(also known as coherence) of matrix [Xi,X 2 ], as well as the 
maximum column correlation of matrices X± and X 2 by 


Ml := 

M 2 ■= 


\{x 1M ,x Ul ’)\ 

max Ti-ii-M-IT’ 

\(X2,h,X 2 , h ')\ 

1 < 1 ™?*< L 2 ||* 2 ,< 2 || 2 • 11 * 2 , i 2 ' II 2 ’ 


uq := max 

1<Zi<Li,1</ 2 <-£/2 


|(*1,<i,*2,< 2 )| 

1 * 1 ,hh ■ \\ X 2,l 2 \ 


(45) 

(46) 

(47) 


taking p = max{^o, A*i, M 2 }. then if 


M < 


1 

Li + L 2 - 1 ’ 


(48) 


we have 

vo\ Ll+L 2 ([U 1 ,U 2 ])>(l-e) L/2 , 


where the matrices U \, U 2 are orthogonal bases for subspaces 
span(Xi) and span(X 2 ), and L = min{Li,L 2 }, 

L 1 L 2 • p 2 

£ - [1 - (Li - l)fj][l - (L 2 - l)g}' 

Proof of Lemma [3J 

Similar to the previous proof, we use the column-normalized 
versions of matrices Xi and X 2 , which are denoted by 
Xi = [* 1 , 1 , ••• ,x ltLl ],X 2 = [* 2 , 1 , ,* 2 ,l 2 ] in the fol¬ 

lowing proof. Therefore we have = *i,! 1 /||*i,t 1 1| 2 , h = 
1, ■ • • , Lx, and x 2M = * 2 ,< 2 /||* 2,< 2 \\ 2 ,l 2 = 1, ■ • • , L 2 ._ It is 
easy to verify that span(Xi) = span(Xi) and span(X 2 ) = 
span(X 2 ). 

According to the relation of volume and principal angles, 
the volume function voli 1 + i 2 ([C/i, U 2 \) satisfies 

L 

VOl Li_+L 2 ac/ 1 ,f/ 2 ])=n sin 0i(span(Xi), span(X 2 )), 
2=1 


where L = min {Lx,L 2 }, and 0 i (span(X 1 ),span(X 2 )) are 
principal angles of subspaces span(Xi) and span(X 2 ). From 
the definition of principal angles, it can be derived that, given 
any principal angle 0j(span(Xi),span(X 2 )), i = !,■■■ ,L, 
there must exist a pair of vectors it,; and v,, which satisfy 

Ui G span(Xi), ||uj|| 2 = 1, (49) 

Vi G span(X 2 ), ||wi|| 2 = 1, (50) 


and 

cos6» i (span(X 1 ),span(X 2 )) = (51) 

Then according to Lemma [2] the condition in < 148 b implies the 
full-rankness of matrix [Xi,X 2 ], which also means that X\ 
and X 2 are full rank; as a result, there must exist a series of 
coefficients aip, • • • , a \and a 2) i, • • • , a 2j l 2 that are not all 
zero, such that 


L 2 

Ui = Oi.ij*!,!]., Vi = a 2: i 2 x 2} i 2 , (52) 

Z 1 —— 1 z 2 =i 


then we have 

L 1 L 2 

\(ui,Vi)\ l a !7i a 2,i 2 | • M 

Zi = l z 2 —1 


< 


\ 


Li ■ ^2 l a L<il 


1=1 


* i 2 - £ l a 2,i 2 1 2 ' Mj (53) 

\ <2=1 


the last inequality is derived from the Cauchy-Schwarz in¬ 
equality. 

Because ||ttj || 2 = 1, we also have 


1 = ||wi|| 2 > 5Z l a L ; i| 2 “ H \ a i,h a i,h'\^ 


h=l 

L\ 


h±h' 


> E Kiil 2 -(ii-i)Z!K'il 2 -^(54) 


< 1=1 


< 1=1 


the last inequality is also derived from the Cauchy-Schwarz 
inequality. Then we can know from (f54t that, 


r. i«i,hi 2 < — 

< 1=1 

similarly we also have 

l 2 

T. \ a 2,h | 2 < 


1 


\-{Lx-l)p' 

1 

1-(L 2 -1)m’ 

combining ( l55l > and ( l56l > with (l53l >. we have 

cos 0;( s P an (-*L)> span(X 2 )) = \(ui,Vi)\ 

< 


(55) 


(56) 


< 2=1 


Li ■ L 2 ■ p 2 


[1 — (Li — l)/i][l — (L 2 — l)/i] ’ 


(57) 


or in another word, 

sin 0 i(span(Xi),span(X 2 )) > 


1 - 


L x ■ L 2 ■ p 2 


[1-(Lx-1)p][1-(L 2 -1)pY 


(58) 


holds for every i, i = 1, • • • , L. If we let e = 

[!_ (Li - L i)^T-(l 2 —ltd ’ Lemma 0 is now proved. 

Proof of Lemma (TJ Theorems [l] and [2} 

First, Lemma Q] is derived based on Lemma [2] Because 
the matrices g[ d ' T ^ in (fTTI) and G 2 in (fl 8 ] > are actually 
composed of different delayed versions of the waveform signal 
g(kT s ), so (| 44 ]> is actually constraining the auto-correlation 
function of g{kT s ). Without loss of generality, we assume 
the signal g(kT s ) to be wide sense stationary, then we can 
approximately have 


\Rs(r ii i i - T it i. 


thus (l45l > and < 146b means 


1 ^( 0 ) I 


i, for li 4 


max \Rg(n,li~ T i,h')\ ^ 

1 <U^U'<Li 


\R a m 

Lx + L 2 — 1 


,i = 1,2, (59) 


because we are deriving a sufficient condition, the condition 
in < IT9t will surely be sufficient to ensure (|59| >. the proof of 
Lemma Q] is completed. 
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As for Theorem Q] the proof is based on Lemma [3] sim- 
ililarly, ( l47l ) is equivalently written as: 


max 

l<h<L 1 ,l<l 2 <L 2 


\R 9 ((t 1m +At-T2,i 2 )T s )\ < 


l^ g ( 0 )l 

Li+La-l’ 


then similarly the condition in (El l is sufficient. 
Theorem |2] is obvious and no need to prove. 


B. Proof of Theorems [?] and |7} 

Before proving Theorems [3] and |4] the following two 
lemmas are firstly given as intermediate results. 

Lemma 4: Consider the linear combinations of several ma¬ 
trices: 


H i = ^2 Qi.ii-ffi.ii = 
Zi=l 

and 

1/2 

H> = ^2 Ot2,l 2 H 2} l 2 = 


y ] Qi.ii ^i,ii >''' j y ' Qi.ii ^i.ii 


(60) 


l 2 


' Q2,i 2 h-2 i 2 ! ■'' j 'y ] Q2 ,i 2 h-2 


(K) 

h 


(61) 

with Hi G C MxK and H 2 =G C MxK , and Wh^h = 

■■■ = ll^iII 2 ’ ll^llb = • • • = ll^lb; we let 


/i 1 = max 


t^S^i iiOiiO= 


p 2 = max 

l<l 2 jtl 2 ’<L 2 




and let 


fiQ = max 

l</l<Li,l</2 <1/2 




l/i 


(fc) 


i,ii 11211 ^ 2 , ; 2 112 


(fe') 


Taking 


B = max{/x 0 ,/xi,/x 2 }, 


(62) 


(63) 


(64) 


(65) 


M?i/iim;/ 1 ii2,--,mi ; /ii^H 2 


.a) 


W/iifcWi 


sumed = ••• = H/i^|| 2 , then the matrices Hi can 

be equivalently written into 


,Wi 


. Because we have as- 


H 1 = 




( 68 ) 


where = a ljh ■ \\h[ k) h \\ 2 , k = l,--- ,iC 
According to the definition of singular value decomposition, 
we have 


min{M,K} 


H 1= £ 




(69) 


Then the left singular vectors of matrix Hi (or H 2 ), cor¬ 
responding to its Ki largest singular values, namely those 
M 11 , • • • , u\ t K-i, can be written as 

= —ffi«i >r , r = 1,• ■ • , Ki, (70) 

0"l,r 

As a matter of fact, these vectors are actually regarded as 
the basis for each receiver’s signal subspaces, and the basis 
matrices are 


Ul = [Ul.l, • • • , Ul : Ki], U 2 =# : |mi, 2 , • • • , (71) 


On the other hand, according to the definition of prin¬ 
cipal angles, take K mm = ruin {A'|, K 2 }, then for every 
principal angle of subspaces span([/i) and span([/ 2 ), i.e., 
0j(span(t/i), span(f/ 2 )), i = 1, • • • ,K m j n , there must exist a 
pair of vectors Mi^ and u 2 ,%, satisfying 

m m G span(C7i), ||'£ti ,*|| 2 = 1, 

« 2 ,» £ span([/ 2 ), ||« 2,*]|2 = 1 , 


such that 


cos 0i(span(f/i), span([/ 2 )) = |(Mi,i, M2,i)|. 

As a matter of fact, these pair of vectors whose angles are 
principal angles, i.e., Mi^ and M 2j i for i = l,--- ,K m ; n , 
are actually the left and right singular vectors of matrix 
U2U 2 - Therefore, they can be orthogonally represented by the 
orthogonal basis matrices Ui and U 2 , which means, there exist 
a series of coefficients q[ x J, ■ ■ • , and q^ 1 }, • • ■ , q 2 K 2 that 
are not all zero, such that 


then we have the same result as in (El l and (Et . 

Lemma 5: Consider the matrices Hi and H 2 in (l60l > and 
(ED. if there exists 1 < Z* < Li and 1 < i| < L 2 , such that 
= H 2i i*\ we let 


Fo 


max 

1 <k,k'<K 

1<Z<^i,Z#ZI,1<Zi<L2,Z 2 #ZS 


(k) 

i,h 


I h 


I h 


W) 
2 ,u 


( 66 ) 


L := min{Li, L 2 }, and let 


Ai = max{/r 0 ,//i,/r 2 } (67) 

then we have the same results as in ( El l and ( [39b . 

Proof of Lemma [4} For simplicity, we take the matrix 
Hi for example. Before we start the proof, for convenience, 
we will use the column-normalized matrices of Hi ^ as our 

main target, namely, we take Hi t i t = [h < 2l l , • • ■ , h'2j] = 


Ki k 2 

Ul,i = ^2q[ r Ju I, r , U2,i = E '? 2 y ? : ) ' u 2,r', 

r— 1 r '=1 

also EE kpil 2 = EEi I92VI 2 = L If we denote 91 4 : = 

[<E, • • • , q[ K i 1 ' > ] T , then we have 

ui,i = Hi ■ - ,v 1}Kl ] ■ diag(cr~E - - - ■ qi,,, 

(72) 

letting Mi ti = [Mr,!,--- ,M liri ] • diag((7 li i,--- ,cr 1 1 Kl ) ■ qi,i, 
then (f72l> are equivalently 

K L\ 

Ui,i = Hi ■vi,i = Y, 4 k i E & LiA k l , (73) 

fc=1 Zi = l 

where v l i = [bj 1 ], • • • , v[ K i ' l } T ■ It is not hard to prove that 
l|ui,i|| 2 = ||diag(cr^i, - - - ,ct 2 Ki ) ■ QiAh < ^ Kl - 


(74) 
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Similarly, if we denote <72,j := [<i \• idiEE’ we have because 
similar results: 


K 


K 


k l 2 

-,(*) 


«2,i = • U2,i = 2 . (75 ) 

fc=l i 2 = l 

ll« 2 ,i|| 2 = l|diag(cr-i, - • - 1 a~ 1 Ki ) ■ q 2> i\\ 2 < ^k 2 (76) 

According to the previous analysis, we have 

= l E E "m^E E test'd 1 < 7? ) 

fc=l fc / = l Z 1 = 1 Z 2 = 1 

< E ■ E • E i"2,i 2 i- E 1 * 4 * E# 1 


- \ & hh\ 2 ■ (84) 


fc=l 


fc=l 


dSTT) becomes 


n ~ 112 

Ul,i 2 


h =1 

therefore we have 

2 


£*$ 


[1 + (Li - 1 )] 71, ( 85 ) 


K 

E ~(fc) 

w i,i 

fc=i 


“ Ef=il«r,/l 2 [l + (ii-lM' 


1=1 

k =1 

l 2 =1 

k' = 1 


K 

< Ei fii 2U 

<=1 \ 

K 

Ewtii 2 

k =1 

l 2 

E ifi2,« 2 u 
<2=1 \ 

K 

K- \ 4 *?\ 2 -» ( 78 ) 

k '=1 

E €’ 

k '=1 


Using the same technique, we can derive 

2 

1 

~~ EEi l« 2 ,( 2 I 2 [1 + (i 2 - l)/r] 


( 86 ) 


(87) 


< E E l«a,i 2 . 

< 1=1 ^,k^2,k 2 


K ■ fi 


(79) |(ui,i,it2,i)| > 


V t 1 + ( Li - ^ I 1 + ( L2 - !)#*] 


-C /2 


holds for every i = l,--- , Amin- The inequality in < 17 8 b 
is based on the Cauchy-Schwarz inequality, and the inequal¬ 
ity in dzu* is based on dH and dm if we take C\ = 

Eq^i l«i ,h\,C 2 = Ef 2 2 =i l« 2 ,z 2 |. Lemma 0] is now proved. 

Proof of Lemma 0 

Similar to the proof of Lem ma l4l we use ( 173b and < 1 77 b to Now ( 1391 ) is proved. We denote 
prove Lemma [5] According to ( 177b . we have 


E E i“2>* 2 1 ■ 

Zl=l,Zl^Z* 12 = 1,12^2 


K 


& 1 ,K\ (J ‘ 2 ,K 2 


■ // ( 88 ) 


i<“r,i>* 2 ,i) 1 = 1E E ^ 4 ^(E«i,<i^. E & 2 ,i 2 h i 2 k l ) 2 )\ 


k=i k'=i 


1=1 


> i«i,i*«2,t*i • E • E i®?* 5 ! - 

k =1 k' = 1 

Li L 2 K K 


E E i fii .L 1 ■ i“ 2 ,*2 i ■ E E 1 * 4 ** 1 • i* 4 k i 5 1 ■ m 

h J 


k =1 k' = 1 


> \ai,qa 2il ,\ ■ I®m I ‘ E 14 VI - 


fc=l fc' = l 




E E Kiil-Kz 2 l- m 

°w*,k 2 


A := — 

sjY.’tU I® 1 ’*! I 2 Ef 2 2 =i |a2,i 2 1 2 

b ■■= E i a i.<iiii^( ) iii 2 • E i“ 2 .' 2 111^2*2 ii 2 

Z 1 =1 Z 2—1 

B* '■= E Kiilll^uJIa ■ E l a 2 ,< 2 lll^ll 2 , 

Zi = l,Zi#Z| l 2 =l,l 2 ^l^ 

Sufficiently, if we let the lower bound of ( 1881 ) be greater than 
the upper bound of (l 79 b . then the right side of d 39 b is less than 
the rig ht side of d 33 b. the corres ponding condition on // is thus 

( 80 ) M < \J ai ' Kl (B+B^-K^ 1) + 3 - \ as in < 5 D- The lemma is 
proved. 


T-, . • . 11 - I, , ■ j j Proof of Theorems l3l and l4l 

Then, the constraint ||m m || 2 = 1 is considered, we can prove xhe result of Lemmali] and Lemma 0 can be directly used 

to derive Theorems 0 and 0 Actually, the matrices H\ and 





h =1 

k= 

2 

1 

K 

2 Li 

K 

2" 

+ 

E &1 ’h 



~ E l“i.iil 2 

E *!5 



h=l 


k = 1 

n=i 

k = l 



H 2 will be replaced by the Hankel matrices x[ Ad T4 and X 2 , 
then the H\^ and Hj i 2 are just S and S 2 ,i 2 \ Therefore 

the expressions (1621164b in Lemma 0 and (l66b in Lemma 0 are 
M (81) actually related with the auto-correlation function of s(k-T s ), 
namely, they can be translated into: 


because the geometrical average is greater than the arithmetic 
average, we have 


l-RsUl,ii + At - r 2 j i 2 + (k — k’)T s 


^ Ef^x lfir.nl 

therefore (18Tb becomes 

ll«i,*lla < E llE a M 1 ®M /l M ) 1 lli+ 


(/o)[At]|. it (k') || 

S l,/l IL ll s 2,Z 2 ll 2 


|fl.(0)| 


(82) 


l( s M 1 [Arl ’ s il ) ' [Ar! >l _ |B.(n,j 1 - r 1M , + (k - k')Ts)\ 


II (fe)[Ar]|| || (fe')[Ar],| 

H S l,il ll 2 ll S l,i 1 ' ll 2 


|fl.(0)| 


I< s 2 fc y4l ) ')l _ I^U 2 ,< 2 -t 2 ,, 2 , + (fc-fe')r B )| 


„( fe ) n„ii=( fe ') 


l«.(o)| 


I-, =1 fc=l 


E Ifir.n I 2 


T,( fc ) 


E ~(l 
V 1 

k =1 


1 The sufficient condition in (182b and (l38b are obtained 

(Li - 1 )fi, (83) through similar techniques in the proofs of Theorem 0 and 

m 















































